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CONVERSION  FACTORS,  U.  S.  CUSTOMARY  TO  METRIC  (SI) 
UNITS  OF  MEASUREMENT 


U.  S.  customary  units  of  measurement  used  in  this  report  can  be  con¬ 
verted  to  metric  (SI)  units  as  follows: 


_ Multiply _ 

feet 

miles  (U,  S.  statute) 


_ B* _ 

0.3048 

1.609344 


To  Obtain 

metres 

kilometres 
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VAHM  -  A  VERTICALLY  AVERAGED  HYDRODYNAMIC  MODEL 
USING  BOUNDARY- FITTED  COORDINATES 

PART  I:  INTRODUCTION 

1.  The  mathematical  modeling  of  the  hydrodynamics  of  a  body  of 
water  plus  the  transport  and  dispersion  of  a  conservative  constituent 
within  that  body  involves  the  solution  of  a  set  of  partial  differential 
equations  expressing  the  conservation  of  mass,  momentum,  and  energy  of 
the  flow  field  along  with  a  transport  equation  for  the  constituent. 

These  equations  involve  derivatives  with  respect  to  time  as  well  as  three 
spatial  dimensions.  However,  a  simplification  that  is  often  made  in 
treating  relatively  shallow  bodies  of  water  that  are  well  mixed  over  the 
depth  is  to  vertically  average  the  three-dimensional  (3D)  equations  to 
yield  a  two-dimensional  (2D)  set  for  nearly  horizontal  flows. 

Numerical  Techniques 

2.  Since  the  governing  equations  are  nonlinear,  analytic  solutions 
in  general  cannot  be  found  and  one  is  forced  to  resort  to  numerical  tech¬ 
niques  to  obtain  solutions.  The  two  most  common  such  techniques  are  the 
finite  difference  method  (FDM)  and  the  finite  element  method  (FEM) . 

There  are,  of  course,  both  advantages  and  disadvantages  to  each  of  these 
approaches . 

3.  Perhaps  the  most  often  quoted  advantage  of  the  finite  element 
method  is  that  with  this  approach  physical  boundaries  coincide  with 
computational  net  points.  Therefore,  the  modeling  of  flow  within  an 
irregular  domain  can  be  more  accurately  handled  than  with  the  normal 
finite  difference  method  where  the  approach  is  to  construct  a  rectangular 
grid  over  the  domain,  which  forces  the  boundaries  to  be  represented  in 

a  "stair  stepped"  fashion.  However,  a  disadvantage  of  finite  element 
methods  is  that  they  involve  dense  matrices  rathei  than  the  sparse  ma¬ 
trices  involved  in  finite  difference  methods.  This  results  in  more 
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computational  time  being  required  in  a  finite  element  model  having 
the  same  number  of  mesh  points  as  a  finite  difference  model.  An  ad¬ 
ditional  disadvantage  is  that  the  finite  element  method  is  more  cumber¬ 
some  to  code  into  a  computer  model  than  the  finite  difference  method. 
This  can  be  a  problem  not  only  during  the  development  of  the  model  but 
can  also  increase  the  level  of  effort  required  during  later  model 
modifications . 


4.  Accepting  that  the  finite  difference  method  possesses  an  ad¬ 
vantage  in  simplicity  and  perhaps  computational  costs,  a  logical  question 
is  whether  or  not  one  can  develop  ways  to  circumvent  the  major  disad¬ 
vantage  of  having  to  represent  irregular  boundaries  in  a  "stair  stepped" 

fashion.  One  such  technique  which  has  been  developed  by  Thompson, 
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et  al.  ’  ’  involves  the  use  of  boundary-fitted  coordinates.  Thompson's 
method  generates  curvilinear  coordinates  as  the  solution  of  two  elliptic 
partial  differential  equations  with  Dirichlet  boundary  conditions,  one 
coordinate  being  specified  to  be  constant  on  the  boundaries,  and  a  dis¬ 
tribution  of  the  other  specified  along  the  boundaries.  However,  the 
numerical  computations  to  solve  the  governing  flow  equations,  as  well  as 
computations  for  the  solution  of  the  coordinate  system,  are  not  made  in 
the  physical  curvilinear  coordinate  system  but  rather  are  made  on  a 
rectangular  grid  with  square  mesh  spacing. 


Purpose  and  Scope 


5.  Since  the  early  to  mid  1960's,  many  finite  difference,  plus  a 
few  finite  element,  computational  models  for  vertically  averaged  flows 
have  been  developed . ^  The  purpose  of  this  report  is  to  describe 
the  development  of  a  new  vertically  averaged  hydrodynamic  model  which  is 
fully  coupled  with  the  water  salinity  through  its  influence  on  the  water 
density.  The  finite  difference  method  of  solution  is  employed  but, 
unlike  the  previously  developed  models,  solutions  are  obtained  on  a 


boundary-fitted  coordinate  system  to  provide  an  accurate  representation 
of  boundary  geometry. 

6.  The  first  part  of  the  report  summarizes  Thompson's  method  for 
computing  boundary-fitted  coordinates.  A  portion  of  this  discussion  has 
been  taken  from  a  previous  Independent  Laboratory  Inhouse  Research  (ILIR) 

g 

report  by  Johnson  and  Thompson.  The  second  part  of  the  report  presents 
the  basic  equations  to  be  solved  and  a  discussion  of  their  transformation 
in  a  fully  conservative  form  from  the  physical  plane  to  a  transformed 
rectangular  plane,  wherein  computations  are  made.  The  third  art  then 
deals  with  the  numerical  aspects  of  the  solution  scheme  and  presents  the 
difference  equations  to  be  solved,  along  with  associated  boundary  con¬ 
ditions.  The  final  part  of  the  report  describes  the  computer  model  as 
it  is  developed  to  date  and  presents  results  from  three  applications 
that  demonstrate,  in  a  qualitative  sense,  that  the  model  is  behaving 
properly. 


PART  II:  ASPECTS  OF  GENERATING  BOUNDARY-FITTED 
COORDINATE  SYSTEMS 


7.  Thompson's  work  on  the  generation  of  boundary-fitted  coordi¬ 
nates  cun  be  found  in  References  1,  2,  and  3.  The  discussion  below  is  a 
summary  of  the  more  important  theoretical  aspects  of  the  subject. 

The  Basic  Idea 

8.  Suppose  one  is  interested  in  solving  a  differential  system 
involving  two  concentric  circles,  such  as  shown  in  Figure  1,  where 

r  =  constant  =  r) ^  on  the  inner  circle  and  r  =  constant  =  n^  on  the 
outer  circle  and  0  varies  monotonically  over  the  same  range  over  both 
the  inner  and  outer  boundaries,  i.e.,  0°  to  360°. 

9.  A  cylindrical  coordinate  system  is  the  obvious  choice  since  a 
coordinate  line,  i.e.,  a  line  of  constant  radius,  coincides  with  each 
boundary.  If  one  now  pulls  the  interior  region  between  the  two  circles 
apart  at  0=0°  (or  0  =  360°)  and  folds  outward,  it  is  easy  to  visu¬ 
alize  the  region  Dj  becoming  the  rectangular  region  D2  . 

10.  The  general  boundary-fitted  system  is  completely  analogous 
to  the  system  discussed  above.  In  Figure  2  the  curvilinear  coordinate, 
ij  ,  is  defined  to  be  constant  on  the  inner  boundary  in  the  same  way  that 
the  curvilinear  coordinate,  r  ,  is  defined  to  be  constant  on  the  inner 
circle  in  the  cylindrical  coordinate  system.  Similarly,  r|  is  defined 
to  be  constant  at  a  different  value  on  the  outer  boundary.  The  other 
curvilinear  coordinate,  £  »  defined  to  vary  monotonically  over  the 
same  range  on  both  the  inner  and  outer  boundaries,  as  the  curvilinear 
coordinate,  0  ,  varies  from  0  to  2n  around  both  the  inner  and  outer 
circles  in  cylindrical  coordinates.  It  would  be  just  as  meaningless  to 
have  a  different  range  for  £  on  the  inner  and  outer  boundaries  as  it 
would  be  to  have  0  increase  by  something  other  than  2n  around  one  of 
the  circles  in  cylindrical  coordinates.  It  is  this  fact  that  £  has  the 
same  range  on  both  boundaries  that  causes  the  transformed  field  to  be 
rectangular.  Note  that  the  actual  values  of  the  coordinates,  r)  and 
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4  ,  are  irrelevant,  in  the  same  way  that  r  and  0  may  be  expressed  in 
different  units  in  cylindrical  coordinates. 

11.  Nov/  that  the  values  of  the  coordinates,  r)  and  4  »  have  been 
completely  specified  on  all  the  boundaries  of  a  closed  field,  it  remains 
to  define  the  values  in  the  interior  of  the  field  in  terms  of  these 
boundary  values.  Such  a  task  immediately  calls  to  mind  elliptic  partial 
differential  equations,  since  the  solution  of  such  an  equation  is  com¬ 
pletely  defined  in  the  interior  of  a  region  by  its  values  on  the  boundary 
of  the  region.  Thus  if  the  coordinates,  4  and  H  *  are  taken  as  the 
solutions  of  any  two  elliptic  partial  differential  equations,  say 
L(4)  =  0  »  D(n)  -  0  ,  where  L  and  D  represent  elliptic  operators, 
tnen  4  and  H  will  be  determined  at  each  point  in  the  interior  of 
the  field  by  the  specified  values  on  the  boundary.  One  condition 
must  be  put  on  the  elliptic  system  chosen  since  the  same  pair  of  values 
(4,n)  must  not  occur  at  more  than  one  point  in  the  field  or  the  co¬ 
ordinate  system  will  be  ambiguous.  This  condition  can  be  met  by 
choosing  elliptic  partial  differential  equations  exhibiting  extremum 
principles  that  preclude  the  occurrence  oc  extrema  in  the  interior  of 
the  field. 


Mathematical  Development 

12.  From  the  discussion  above,  a  logical  choice  of  the  elliptic 
generating  system  is  Poisson's  equation.  Thus,  based  upon  Figure  2,  the 
basic  problem  is  to  solve 


4  +4 

xx  yy 


=  P 


0  + 
xx 


n  =  Q 
yy 


(i) 


with  boundary  conditions, 
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(2) 


i  =  £j(xky)  on  Tj 
H  =  constant  =  r)j  °n  Tj 

£  ~  42(x,y)  on  r2 

r)  =  constant  =  r)2  on  ^2 

The  arbitrary  carve  joining  and  P^  in  the  physical  plane  specifies 

a  branch  cut  for  the  multiple-valued  function,  £(x,y)  •  Thus  the  values 
of  the  coordinate  functions  x(£,r|)  and  y(£,n)  are  equal  along  P^ 
and  P.  ,  and  these  functions  and  their  derivatives  are  continuous  from 
r3  to  P^  .  Therefore  boundary  conditions  are  neither  required  nor 

allowed  in  P_  and  T,  . 

3  4 

13.  The  functions  P  and  Q  may  be  chosen  to  cause  the  coordi¬ 
nate  lines  to  concentrate  as  desired.  As  discussed  in  Reference  1, 
negative  values  of  Q  result  in  a  superharmonic  solution  and  cause  r) 
lines  to  move  toward  the  r)"liIie  having  the  lowest  vaue  of  I]  ,  while 
positive  values  have  the  opposite  effect.  Considering  the  £  solution 
to  be  superharmonic  results  in  the  interior  of  the  £  =  constant  lines 
being  rotated  in  a  clockwise  direction  in  the  physical  plane,  whereas  if 
the  £  equation  is  subharmonic,  i.e.,  P  is  positive,  the  lines  are 
rotated  in  the  counterclockwise  direction. 

2 

14.  The  form  of  these  functions  incorporated  by  Thompson,  based 
upon  much  computer  experimentation,  is  that  of  decaying  exponentials. 

For  example,  let  Q  be  taken  as 


Q  = 


exp  (-  d  n  -  ni  ) 


where  a  and  d  are  constants,  and  r)^  is  some  specified  q-line. 
This  function  reaches  its  maximum  magnitude  on  the  r)^  line  and  decays 
away  from  that  line  on  either  side  at  a  rate  controlled  by  d  . 
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15.  This  function  would  cause  r|— 1  ines  to  concentrate  on  one  side 
of  the  r)^-.'ine  and  to  move  away  from  the  other  side.  If,  however,  a 
sign-changing  function  is  incorporated  so  that 

Q  =  -  a  sgn  (n  -  ni)  exp  (-  d|n  -  ni| ) 

where  sgn(x)  is  simply  the  sign  of  x  ,  the  rj  —  1  ines  will  concentrate 
on  both  sides  of  the  r;^-line.  In  a  similar  fashion,  it  is  possible  to 
cause  concentration  of  rj ^  —  1 ine s  near  a  point  (£^,0^)  with  the  function 

Q  =  -  a  sgn  (n  -  Hi)  exp  £  -  d^(£-£i)2  +  (rl-ni)2  J 

Finally,  concentration  near  more  than  one  line  and/or  point  is  achieved 
by  writing  Q  as  a  sum  of  functions  of  the  above  form.  In  this  case 
the  attraction  amplitude  a  and  the  decay  factor  d  may  be  different 
for  each  line  or  point  of  attraction.  The  decay  factor  should  be  large 
enough  to  cause  the  effects  of  each  attraction  line  or  point  to  be  con¬ 
fined  essentially  to  its  immediate  vicinity.  Thompson  has  found  that 
attraction  amplitudes  of  100  are  moderate,  10  is  weak  and  1000  is  fairly 
strong.  A  decay  factor  of  1.0  causes  the  effects  to  be  confined  to 
a  few  lines  near  the  attraction  source,  while  0.1  gives  a  fairly  wide¬ 
spread  effect.  Control  of  £-lines  is  accomplished  by  an  analogous 
form  of  the  function  P  .  Such  control  is  useful  to  improve  grid  spacing 
and  configuration  when  complicated  geometries  are  involved. 

16.  Since  all  numerical  computations  are  to  be  performed  in  the 
rectangular  transformed  plane,  it  is  necessary  to  interchange  the  de¬ 
pendent  and  independent  variables  in  Equation  1.  Using  the  relations 


-  yj 

=  -yj 

nx  =  -yi/J 


ny  =  xc/j 
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'  i  ii  — aitfii— ‘iiiim1- 


«xx  =  <{xyin  *  Vnn)/J  ‘  (&t  +  «x'>xJn)/J 

4yy  =  -  (Vnn  +  VW/J  ‘  (WJn  +  $tw 

nxx  =  -  («xy5i  +  Vjn)/J  '  (t*Vi  *  nxJn)/J 

r)  =  (n  x.  +  £  xtt)/J  -  (£  n  J«.  +  n2J  )/J 

yy  y -on  y  u  y  y£  y  n 

equation  1  becomes 

“x«  '  2Px{n  ‘  Yxnn  +  j2(Px{  +  Qxn)  *  0 

(3a) 

“y«  ■  2PySn  +  vynn  *  j2(Py£  +  Qyn)  =  0 

wh;re 


2  ,  2 
a  =  x  +  y 

n  yn 


P  =  »ixn 


yey 


4yn 


+ 


2 

y4 


J  =  Jacobian  of  the  transformation  =  x 


x«y 


ny£ 


with  the  transformed  boundary  conditions 


(3b) 


11 


x  =  f  1(4»n1'> 


on 


y  =  g1(t.n1)  on  r* 

(4) 

x  =  f2^,lV  °n  P2 

y  =  82(4,n2)  on  r2 

Again  considering  Figure  2,  the  functions  fj(|,rij)  ,  » 

f2(C,n2)  »  anfl  g2^,r^  are  sPec^^e<^  by  the  known  shape  of  the  con¬ 
tours  Tj  and  and  the  specified  distribution  of  4  thereon.  Al¬ 

though  the  new  system  of  equations  is  more  complex  than  the  original 
system,  the  boundary  conditions  are  specified  on  straight  boundaries  and 
the  coordinate  spacing  in  the  transformed  plane  is  uniform.  Computa¬ 
tionally,  these  advantages  far  outweigh  any  disadvantages  resulting  from 
the  extra  complexity  of  the  equations  to  be  solved. 

17.  The  boundary-fitted  coordinate  system  so  generated  has  a 
constant  q~line  coincident  with  each  boundary  in  the  physical  plane. 

The  4-lines  may  be  spaced  in  any  manner  desired  around  the  boundaries 
by  specification  of  x,y  at  the  equispaced  4“P°i.nts  on  the  T*  and 
F*  lines  of  the  transformed  plane. 

18.  The  rectangular  transformed  grid  is  set  up  to  be  the  size 
desired  for  a  particular  problem.  Since  the  values  of  4  and  q  are 
meaningless  in  the  transformed  plane,  the  q  lines  are  assumed  to  run 
from  1  to  the  number  of  q  lines  desired  in  the  physical  plane.  Like¬ 
wise,  the  4  lines  are  numbered  1  to  the  number  specified  on  the  bound¬ 
aries  of  the  physical  plane.  The  grid  spacing  in  both  the  4  an<^  H 
directions  of  the  transformed  plane  is  taken  as  unity.  Second  order 

central  difference  expressions  are  used  in  Thompson's  coordinate  genera- 
2 

tion  code,  TOMCAT,  to  approximate  all  derivatives  in  Equations  3a  and 
3b.  The  resulting  set  of  nonlinear  difference  equations,  two  for  each 
point,  are  solved  in  TOMCAT  by  accelerated  Gauss-Seidel  (SOR)  iteration 
using  overrelaxation.  Some  discussion  of  this  technique  is  presented  in 
Reference  2. 
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19.  The  same  procedure  may  be  extended  to  regions  that  are  more 
than  doubly  connected,  i.e.  have  more  than  two  closed  boundaries,  or 
equivalently,  more  than  one  body  within  a  single  outer  body.  A  river 
reach  containing  more  than  one  island  would  be  an  example. 

Types  of  Boundary-Fitted  Coordinate  Systems 

20.  Previous  discussion  of  the  generation  01  boundary-fitted 
coordinates  has  centered  around  the  idea  of  using  branch  cuts  to  reduce 
multiply  connected  regions  to  simply  connected  ones  in  the  transformed 
plane.  Thompson's  TOMCAT  code  employs  such  branch  cuts.  The  other  type 
of  coordinate  system  transformation  available  leaves  the  multiplicity  of 
the  region  unchanged.  In  this  case,  bodies  in  the  interior  of  the 
physical  field  are  transformed  to  rectangular  slabs  or  even  slits  in  the 
transformed  plane.  In  the  case  of  slits,  the  physical  coordinates  and 
solution  variables  generally  have  different  values  at  points  on  the  two 
sides  of  the  slit,  even  though  such  points  are  coincident  in  the  trans¬ 
formed  plane.  This  does  not  introduce  any  approximations,  but  simply 
adds  a  little  more  bookkeeping  to  the  code.  Fields  with  more  than  one 
body  in  the  interior  simply  result  in  a  like  number  of  slabs  and/or 
slits  in  the  transformed  plane. 

21.  Different  types  of  transformation  may  be  more  appropriate  for 
different  physical  configurations.  Generally,  the  slit/slab  form  is 
more  appropriate  for  channel-like  physical  configurations  having  bodies 
in  the  interior,  while  the  branch  cut  form  works  particularly  well  for 
"unbounded"  regions  involving  external  flow  about  bodies  and  for  regions 
having  an  outer  boundary  that  forms  a  continuous  circuit  without  pro¬ 
nounced  corners  around  the  field.  The  slab  is  generally  superior  to  the 
slit  unless  the  boundary  has  a  sharp  point.  The  case  of  a  single  channel 
without  any  interior  bodies  would  be  the  same  in  either  form. 

Data  Required  for  Generation  of  Boundary-Fitted  Coordinates 

22.  The  basic  input  or  data  required  to  generate  a  boundary-fitted 
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coordinate  system  are  the  physical  coordinates  of  points  on  the  bound¬ 
aries.  This  will  be  discussed  in  more  detail  in  PART  V  in  connection 
with  the  applications  presented. 


Computer  Time  Required  for  Generation  of 
Boundary-Fitted  Coordinates 

23.  The  computing  cost  for  generating  a  boundary-fitted  coordinate 
system  is  trivial.  Approximately  3  sec  of  CPU  time  on  a  CRAY  I  computer 
were  required  to  generate  the  coordinate  system  shown  in  Figure  3.  It 
might  be  noted  that  no  coordinate  control  was  employed.  The  use  of  such 
control  would  result  in  a  slight  increase  in  computational  time. 
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PART  III:  BASIC  HYDRODYNAMIC  EQUATIONS 


24.  The  Navler  Stokes  equations  express  the  conservation  of  mass 
and  momentum  of  a  flow  field  and  are  the  basic  governing  equations  for 
the  solution  of  any  fluid  dynamics  problem.  Written  in  tensor  notation 
these  equations  are 


Continuity: 


3pu. 

3x . 

1 


0 


(5) 


-3P 


3pui  3(pu.u.) 

Momentum:  -v.-—  +  — s „ 

-  3t  3x.  3x. 

J  i 


Pg.  -  2t.  ,kn  Puk  ♦ 


3T. 

1 


3x. . 

J 


(6) 


where 


e .  . .  = 
ijk 

fi.  = 
J 

T.  .  = 
ij 

M  = 

6.  .  = 
ij 


and  where 


water  density 
time 

tensor  notation  for  velocity 

tensor  notation  of  spatial  coordinate 

acceleration  of  gravity 

cyclic  tensor 

Coriolis  parameter 

laminar  stress  tensor 

molecular  eddy  viscosity 

Kronecker  delta 


(7) 


represents  the  viscous  molecular  stress  arising  as  a  result  of  the  con 
tinuum  approach.  All  symbols  used  are  defined  in  Appendix  A.  It  will 
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fc.'  *&lw 


be  recalled  from  tensor  theory  that  repeeved  indicies  imply  a  summation 
and  also  that  c. . u  in  the  Coriolis  term  is  the  cyclic  tensor  defined 

1  -  for  an  even  permutation  of  ijk 
1  -  for  an  odd  permutation  of  ijk 
0  -  otherwise 

In  addition,  the  Kronecker  Delta,  6ij  ,  is  defined  as 

1  -  if  i  =  j 

6 . .  =0  -  otherwise 
ij 

25.  In  addition  to  the  above  equations,  a  conservation  of  mass 
equation  must  also  be  written  for  any  constituent  being  transported. 
Such  an  equation  for  the  salinity  becomes 


Salinity: 


(8) 


This  equation  states  that  the  salinity  can  change  as  a  result  of  advec- 
tion  by  the  flow  field  and  molecular  diffusion. 

26.  Since  the  salinity  is  coupled  to  the  flow  equations  through 
its  influence  on  the  density,  one  additional  equation  remains  to  be 
written  in  order  to  close  the  system.  An  equation  of  state  expressing 
the  density  as  a  function  of  the  temperature  and  salinity  must  be 
employed . 


Equation  of  State:  p  =  p(T,s)  (9) 

With  the  closure  of  the  system,  there  exists  six  equations  to  be  solved 
for  the  six  unknowns*;  density  -p  ,  three  velocity  components  -u,v,w  , 
pressure  -p  ,  and  salinity  -s. 
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Time  Averaging  for  Turbulent  Flows 


27.  The  above  equations  written  with  molecular  values  of  viscosity 
and  diffusivity  are  only  applicable  in  a  practical  sense  to  laminar  flow 
fields  where  the  flow  does  not  exhibit  random  irregular  fluctuations  in 
time.  However,  most  fluids  in  motion  exhibit  such  fluctuations  and  are 
referred  to  as  turbulent  flows. 

28.  following  Reynolds,  the  approach  normally  taken  to  make  the 
equations  applicable  to  turbulent  flows  is  to  assume  that  the  dependent 
variables  are  composed  of  an  average  time-varying  component  plus  a  small 
randomly  varying  component  about  the  average  value.  This  is  illustrated 
below. 


Thus,  one  writes 


u.(x,y,z,t)  =  u,(x,y,z,t)  +  u^(x.y,z,t) 


where 


t+At/2 


u . 
r 


1_ 

At 


/ 


(x ,y , z , t)  dt 


t-At/2 


A 


:[ 
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and 


kt 

t+At/2 

u^(x,y,z,t)  dt  =  0 

t-At/2 

u!  -  deviation  between  instantaneous  velocity  and  time-averaged 
velocity 

u^  =  time-averaged  velocity 
At  =  time  step 

With  all  the  dependent  variables  written  in  the  form  above,  substitution 
into  Equations  5,  6,  and  8  and  then  integration  over  the  time  increment 
At  produces  the  same  form  of  the  previous  equations,  but  now  written 
with  the  time-averaged  C'\.ifouein.s  as  the  dependent  variables,  plus  the 
additional  terms 

t+At/2 

~T  f  u’.u’.  dt 
At  J  i  J 

t-At/2 

and 

t+At/2 

h  f  s'uidt 

t-At/2 

where  s’  =  deviation  between  instantaneous  and  time-averaged  salinity. 

29.  The  first  term  is  referred  to  as  the  turbulent  Reynolds  stress, 
since  the  high  frequency  turbulent  fluctuations  manifest  themselves  as 
viscous  stresses  acting  on  the  average  component  of  flow.  Using 
Boussinesq's  concept  of  eddy  viscosity,  the  first  term  is  written  as 
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t+At/2 


1_ 

At 


/ 

t-At/2 


u'.u'. 
i  J 


dt 


summation  over 


i) 


In  analogy  with  the  laminar  flow  case,  £..  is  referred  to  as  the 

ij 

turbulent  or  eddy  viscosity  tensor. 

30.  In  a  similar  fashion,  the  second  term  above,  which  arises 
from  the  time  averaging  of  the  salinity  equation,  is  commonly  written  as 


t+At/2 

it  /  s'ui dt  =  Au  ti 

t-At/2 


where  A_  is  called  the  "eddy  diffusivity  tensor"  and  s  is  the  time- 
averaged  salinity. 

31.  The  equations  commonly  applied  to  turbulent  flow  problems  can 
now  be  written  as 


Cq^ipuity:  §*:  +  =  0 


(10) 


Momentum: 


3pu.  3(pu  u  ) 


at 


9x . 
J 


ap 

aT  + 

1 


2e.  ..n.pu.  +  „ 

ljk  y  k  9x, 


jl  iA9xj  9xi/_ 
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Salinity: 


(12) 


9s  ^SUi  _  9  / .  9s  \ 

9t  9x^  -  9x^  ij  9Xj  J 

Equation  of  State:  p  =  p(T,s)  (13) 


where 

p  =  time-averaged  water  density 
P  =  time-averaged  pressure 

and  where  the  assumption  has  been  made  that  the  eddy  coefficients  are 
much  larger  than  the  molecular  values;  i.e., 

e . .  >>  p 

ij 


Depth  Averaging  for  Nearly  Horizontal  Flow 

32.  A  solution  of  the  above  set  of  equations  constitutes  a  fully 
time  varying,  three-dimensional  model  of  the  flow  and  salinity  fields. 
However,  when  modeling  nearly  horizontal  flow  in  relatively  shallow  and 
well-mixed  water  bodies  the  usual  approach  is  to  employ  a  spatial  averag¬ 
ing  to  yield  a  two-dimensional  model. 

33.  The  basic  assumption  in  the  spatial  averaging  of  the  three- 
dimensional  equations  is  that  the  dependent  variables  can  be  represented 
by  an  average  value  over  one  or  more  of  the  spa  ial  coordinates  plus 
some  small  random  deviation;  e.g.,  the  velocity  would  be  written  as 


where 


u.  -  u.  +  u', 
ill 


(14) 
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x.+Ax./2 


1  1  1 
—  f  u.  dx. 

*i  J 


x.-Ax./2 
1  l 


x.+Ax./2 
1  1 


—  f  u!  dx.  =  0 

x  J  11 

*•  __  A /  A 


x.-Ax./2 
i  l 


=  time-  and  space-averaged  velocity 
Ax^  =  spatial  step 

u'  =  deviation  between,  time-averaged  velocity  and  time-  and  space- 
averaged  velocity 

In  an  x  ,  y  ,  z  coordinate  system  (with  x  referring  to  the  longi¬ 
tudinal;  y  ,  the  lateral;  and  z  ,  the  vertical),  if  i  =  2  ,  the  inte¬ 
gration  is  over  the  width  and  a  width-averaged  model  results.  However, 
if  i  =  3  ,  the  integration  is  taken  over  the  depth  and  a  depth-averaged 
model  will  result.  Many  depth-averaged  models  nave  been  developed  since 
Leendertse ' s^  work,  whereas  laterally  averaged  models  have  only  been 
developed  over  the  past  five  years  or  so.  If  the  integration  is  per¬ 
formed  over  the  complete  cross  section,  a  one-dimensional  model  with 
variations  allowed  only  in  the  longitudinal  direction  results. 

34.  As  was  done  in  the  time-averaging  of  the  instantaneous  equa¬ 
tions,  expressions  such  as  Equation  14  are  substituted  into  the  turbulent 
time-averaged  equations  to  yield  a  set  of  equations  with  the  time- 
averaged  and  spatially  averaged  components  of  the  flow  and  salinity  as 
dependent  variables  plus  the  additional  terms 


x.+Ax./2 

l  i 


r- —  I  U  U  dx  . 

ix.  J  1  J  1 


x.-Ax./2 

i  l 
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and 


x . +Ax . /2 
1  i 


x . -Ax . /2 

1  l 


dx . 

l 


As  in  the  time-averaging  case,  these  terms  are  normally  approximated  by 


1 

Ax . 


x.+Ax./2 
i  „  i 


/ 


-Ax./2 


u !  u dx . 
i  J  i 


e!  . 

ij 


+ 


) 


and 


x .+Ax . /2 

i  1  r1 

t —  /  s  ’  u !  dx .  =  A 

Ax.  I  ii 

1  J 

x.-Ax./2 

i  l 


X 


where  e! .  and  A! .  are  referred  to  as  "eddy  dispersion  coefficients"  by 
g  iJ  iJ 

Holley  to  distinguish  them  from  the  turbulent  eddy  diffusion  coeffi- 

**s/ 

cients  arising  from  the  time  averaging,  and  s  is  the  time-averaged  and 
and  spatially  averaged  salinity. 

35.  The  resulting  spatially  averaged  equations  take  different 
forms,  depending  upon  whether  the  averaging  is  performed  over  the  depth 
or  the  width.  For  the  depth  averaged  case,  the  equations  below  are 
obtained.  It  should  be  noted  that  the  Boussinesq  approximation  has  been 
made  which  removes  the  effect  of  density  variations  in  all  terms  except 
those  multiplied  by  the  acceleration  of  gravity. 

„  Si)),  3(uh)  ,  3(vh) 

Continuity:  £  +  — =  0 
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x-momentum: 


3(hu)  +  3(hu2)  +  3(huv)  h  3P 

3t  3x  3y  p  8x 


3  hD 


xx 


3u  \ 
3x  / 


3x 


3  hD 


_X£ 


3y 


3u  \ 
3y  / 


+ 


XB  +  fhv 
X 


(16) 


y-momentum: 


3(hv)  +  3(huv)  3  (hv"1) 
3t  3x  3y 


h  3P 
P  3y 


■K,  1 ) 

3x 


I  (hD. 


3v  \ 

yy  3y  / 


3y 


+ 


fhu 


(17) 


Salinity: 


Ws)  +  3hus  +  30WJ)  _  3K  H)  ,  itsM) 

3t  3x  3y  3x  3y 


(18) 


The  equation  of  state  relating  the  water  density  to  the  salinity  and 
water  temperature  (assumed  constant)  has  been  taken  from  Leendertse* 
and  is  given  as 


p(s,T)  =  1000  +  ALO  *  P0)  (19) 

where 

AL  =  1779.5  +  11.25T  -  0.0745T2  -  (3.80  +  0.01T)s 
ALO  =  0.6980 

PO  =  5890.0  +  38T  -  0.375T2  +  3s 
36.  In  the  above  equations  the  surface  wind  shear  is 
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(20) 


T 

s 

X 


c  2 

—  p  v  cos  a 
paw 
o 


and  the  bottom  shear  is 


c  2 

—  p  v  sm  a 
p  aw 
o 


The  coriolis  parameter,  f  ,  is  computed  from 

f  =  2ui  sin  A 
e 


(21) 


(22) 


(23) 


(24) 


where  =  earth's  angular  velocity  and  A  is  the  angle  of  latitude  of 
the  center  of  the  area  being  modeled. 

37.  In  order  to  finalize  the  above  system  of  equations  it  remains 
to  couple  the  salinity  computations  with  those  of  the  flow  field.  This 
is  accomplished  in  the  following  manner.  Assuming  that  the  pressure  is 
hydrostatic , 


9P 

3z 


"PS 


one  can  determine  the  pressure  at  any  depth  z  from 


pgdz 
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where  the  coordinate  system  is 


Integrating  the  above  eq  one  obtains 


♦ 


P 


pgdz 


z 


where  P  is  the  atmospheric  pressure.  Differentiating  with  respect  to 
the  x-coordinate  yields 


8P_?j  8_  j 

3x  "  ax  ax  j  pgdz 

z 

As  was  done  in  the  continuity  and  momentum  equations  above,  one  now 
assumes  that  the  pressure  and  density  are  composed  of  a  depth  averaged 
plus  a  fluctuating  component.  The  resulting  equation  is  then  integrated 
over  the  depth  to  yield 


h 


3P 

3x 


ap 

h  ^+hgP 


9<t> 

ax 


+ 


1 

2 


h2g 

n  8  ax 


(25) 
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Similarly s 


38.  As  previously  discussed,  the  above  set  of  equations  must  now 

be  transformed  into  a  (4,0)  boundary-fitted  coordinate  system  such  that 

(4,0)  are  the  independent  variables.  The  resulting  set  of  equations  will 

then  be  solved  in  a  transformed  rectangular  plane  as  discussed  in  PART  II. 

In  order  to  accomplish  the  transformation,  the  following  expressions  de- 
12 

rived  by  Thompson  are  utilized. 


[(fVt  ■  CfVn] 

(32) 

[  -  (fx„)t  +  <fVn] 

(33) 

It  should  be  noted  that  these  expressions  are  written  in  a  fully  conserva¬ 
tive  form  which  should  result  in  a  more  accurate  solution  in  highly  ir¬ 
regular  coordinate  systems. 

39.  Using  the  above  expressions  and  assuming  that  the  coordinate 
system  is  time  invariant,  one  can  transform  Equations  27-31  into  the  set 
below. 


27 


Transformed  Equations 


Continuity: 


x-Momentum: 


+ 


y-Mo men turn: 


If  +  T  [  ( uhyn  '  vhxn)f  +  (vhx?  ~  1  -  0 


(34 


9(hu)  ,  _1 

at  j 


( 


hu  y^  -  huvx 


n)?  +  (1’uvx5 '  hu2yE)n  ] 


'  '  K  [  (P‘H  'M,  ]'k[  K)s  -  (*y0n] 

■^[K)t^)J4j(¥4K)t 

-(^[-(%)5  +  K)n]4  +  (>[  -K)f 

K) 


XK  L  {  pavw  cos  “  -  k2  +  fhv  (35) 


8(M  ,  I 


at 


[  (hv%  '  huv^)  +(huvyn  -  hv2xny  J 

■ '^[(pA)n-(‘v‘4]  ~*t  L  "(**"), 

+K)J-I^[  -K\+K)J 


40.  The  above  set  of  equations  constitute  the  set  for  which  a 
numerical  solution  is  sought  on  a  rectangular  grid  with  square  grid 
spacing  (e.g.,  A6  ■  At;  ■  1.0)  .  It  remains,  of  course,  to  specify  proper 
boundary  conditions  along  the  sides  of  the  rectangular  grid.  It  is 
obvious  that  the  transformed  equations  are  more  complicated  than  the 
original  cartesian  forms;  however,  the  advantage  of  being  able  to  make 
computations  on  a  rectangular  grid  far  outweighs  any  disadvantage  result¬ 
ing  from  the  more  complicated  set  of  equations. 
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PART  IV:  NUMERICAL  ASPECTS 


41.  In  order  to  obtain  a  solution  of  the  governing  set  of  Equa¬ 
tions  34-38,  the  method  of  finite  differences  is  employed.  There  are 
many  different  types  of  finite  difference  schemes  that  have  been  employed 
in  numerical  solutions  of  partial  differential  equations.  These  schemes 
range  from  fully  explicit  to  fully  implicit,  with  a  combination  of  an 
explicit-implicit  scheme  being  employed  in  some  cases,  e.g.,  Edinger  and 
Buchak.1  A  similar  scheme  is  employed  here.  Basically,  the  computa¬ 
tional  cycle  will  consist  of  the  following  steps: 

a.  Solve  for  the  water  surface  from  the  continuity  equation 
in  a  fully  implicit  fashion. 

ti.  Using  the  most  recent  values  of  the  water  surface  eleva¬ 
tions,  solve  for  the  u  and  v  velocity  components  from 
the  x  and  y  momentum  equations  in  an  explicit  fashion. 

£•  Solve  for  the  salinity  from  the  salt  transport  equation 
in  an  explicit  fashion. 

d.  Compute  the  density  from  the  equation  of  state,  using  the 
most  recently  computed  salinity  field. 

e.  Step  forward  in  time  and  repeat,  the  sequence. 

Such  a  scheme  as  outlined  above  will  have  the  stability  criterion  asso¬ 
ciated  with  the  speed  of  a  free  surface  gravity  wave  removed;  although, 
diffusive  criteria  as  well  as  the  Torrence  condition  associated  with  the 
speed  of  a  water  particle  remain.  However,  these  criteria  are  not 
normally  overly  restrictive. 


Computational  Grid 

42.  The  grid  upon  which  Equations  34-38  are  solved  is  rectangular 
with  a  grid  spacing  of  A£  =*  An  ®  1  .  The  u  and  v  velocity  compo¬ 
nents  are  computed  at  the  corners  of  each  cell  with  the  water  surface 


elevation,  salinity,  and  density  computed  at  the  center  of  a  cell.  Such 
a  grid  is  illustrated  below. 


The  (x,y)  coordinates  are  specified  at  the  corners,  the  center,  and  also 
at  the  midpoint  of  each  aide  of  a  cell. 

43.  One  might  think  of  the  above  grid  as  a  global  grid.  A  local 
grid  consisting  of  25  points  surrounding  the  (£,n)  point  at  which  compu¬ 
tations  are  being  made  is  utilized  in  writing  the  difference  form  of  the 
governing  equations.  This  grid  is  as  shown  below  when  velocity  computa¬ 
tions  are  being  made. 


With  the  above  grid,  when  u  and  v  are  being  computed  at  the  point 
labeled  c  ,  velocities  are  defined  at  points  NWNW,  NN,  WW,  NENE,  EE, 
SESE,  SS,  and  SWSW,  whereas  the  water  surface  and  salinity  are  defined 
at  NW,  NE,  SE,  and  SW.  When  a  value  of  one  of  the  dependent  variables 
is  needed  at  some  point  where  the  variable  is  not  defined,  an  averaging 
is  performed,  e.g.,  u(W)  =  [u(c)  +  u(WW)]  /  2  . 

44.  Points  in  the  local  grid  are  assigned  to  their  location  in 
the  global  grid  through  functions  of  the  form 

IFCOR(L)  -  I  +  INFCOR(L) 

JFCOR(L)  =  J  +  JNFCOR(L) 

where  the  25  values  of  INFCOR(L)  are 

-1  0  0  0  1 

0-1000 
-1  0  0  0  1 

0-1000 
-1  0  0  0  1 

and  the  values  of  JNFCOR(L)  are 

-1  0-1  0  -1 
0-1  0-1  0 
0  0  0  0  0 

0  0  0  0  0 

+1  0+1  0  +1 

Thus,  as  an  example,  if  one  considers  the  computation  for  u  at  (5,5) 
then  the  value  of  u(WW)  in  the  local  grid  should  correspond  to  the 
value  at  (4,5)  in  the  global  grid.  Using  the  expressions  above,  with 
WW  *  11  from  the  local  grid,  one  obtains 
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46.  As  previously  noted,  the  water  surface  elevations  are  to  be 
computed  using  an  implicit  scheme.  Thus,  in  writing  the  difference  form 
of  the  continuity  equation  all  spatial  derivatives  are  taken  at  the  new 
time  level  (n+1)  .  Equation  34  becomes 


47.  In  the  x  and  y  momentum  equations,  all  terms  are  taken  at 
the  old  time  step  except  the  water  surface  slope  term  which  is  computed 
at  the  new  time  step.  Therefore,  the  difference  form  of  the  x  and  y 
momentum  equations  becomes 


+  Fn  (40') 

c 


c 


(41) 


where  F  and  G  contain  all  other  terms  in  Equations  35  and  36,  respec 
tively.  The  difference  forms  of  F  and  G  are 
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and 
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48.  Consider  the  cell  below 


From  an  inspection  of  Equation  39,  it  can  be  seen  that  one  needs  (uh) 
and  (vh)  on  the  cell  faces  at  time  level  (n+1)  in  order  to  solve  for 
the  water  surface  elevation  at  the  center  of  the  cell.  From  Equations  40 
and  41,  one  can  determine  (hu)n+^  and  (hv)n+^  at  the  cell  comers. 
From  these  values  one  can  determine  values  on  the  faces  by  averaging, 

e.g. . 


Now  if  one  substitues  into  Equation  39  for  the  values  of  (uh)n+*  and 

n_|_l 

(vh)  on  the  faces  (from  Equations  40  and  41  with  appropriate  averag¬ 
ing)  an  equation  containing  only  (j>  at  the  (n+1)  time  level  results. 
This  equation  is  then  solved  fcr  by  using  the  Accelerated  Gauss- 

Seidel  solution  technique. 

49.  After  the  water  surface  elevation  at  the  center  of  each  cell 
is  determined  at  the  (n+1)  time  level,  values  of  un+^  and  v°+*  at 
the  cell  corners  are  explicitly  determined  from  Equations  40  and  41 
using  the  new  <j>’s  at  the  (n+1)  time  level.  It  might  be  noted  that 
the  expressions  for  F  and  G  in  Equations  42  and  43  are  only  computed 
once  during  each  time  step.  These  values  are  then  used  in  first  the 
iteration  on  the  water  surface  and  then  in  the  velocity  computations. 

50.  In  the  computation  of  <j>  ,  u  ,  and  v  ,  the  density  is 
taken  at  the  old  time  level.  Its  value  at  the  new  time  level  is  computed 
from  the  equation  of  state  relating  the  density  to  the  salinity  at  the 
new  time  level.  New  salinities  are  computed  from  an  explicit  representa¬ 
tion  of  the  salt  transport  Equation  37.  The  difference  form  becomes 
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Step  1:  Compute  the  terms  labeled  F  and  G  in  t''e  x  and 
y  momentum  equations  at  the  cell  corners  fi^.j  Equa¬ 
tions  41  and  43. 

Step  2:  Using  the  Accelerated  Gauss-Seidel  solution  technique, 
implicitly  solve  for  the  water  surface  elevation  at 
the  center  of  each  cell. 
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Step  3:  Using  the  new  water  surface  elevations  and  values 

of  F  and  G  from  Step  1,  explicitly  solve  for  the 
velocity  components,  u  and  v  ,  from  Equations  AO 
and  Al. 

Step  A:  Explicitly  compute  the  salinity  field  from  Equation  AA. 

Step  5:  Using  the  new  salinities  from  Step  A,  compute  the 
water  density  from  the  equation  of  state. 

Step  6:  Update  all  arrays,  increment  the  time,  and  return  to 
Step  1. 


52.  This  solution  scheme  removes  the  gravity  wave  stability 
criterion  from  consideration,  although  it  should  be  noted  that  other 
stability  criteria  still  control  the  size  of  the  computational  time  step 
allowed.  These  criteria  have  not  been  derived  for  the  transformed  equa¬ 
tions;  however,  for  the  cartesian  form  of  the  equations  they  are 
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In  other  words,  the  time  step  must  be  small  enough  so  that  a  fluid 
particle  does  not  move  more  than  one  grid  spacing  during  the  time  step. 
This  basic  criterion  is  not  nearly  as  severe  as  the  gravity  wave  criterion 
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for  most  practical  problems. 


Boundary  Conditions 


53.  Three  types  of  boundaries  are  allowed  in  VAHM;  walls,  oceans, 
and  rivers.  Wall  boundaries  are  characterized  by  the  specification  of  a 
no-slip  condition,  i.e.,  the  velocity  components  u  and  v  are  set  to 


be  zero  at  walls.  Although,  physically,  the  flow  must  be  zero  at  a 
solid  boundary,  slip  conditions  on  the  velocity  at  a  wall  often  give 
more  realistic  results  if  the  grid  spacing  is  too  large  near  the  wall. 
Slip  conditons  would  be  implemented  by  setting  the  normal  component  of 
the  velocity  equal  to  zero  with  the  tangential  component  computed  from 
the  expression  for  zero  vorticity.  At  the  present  time,  only  the  no¬ 
slip  condition  is  allowed  in  VAHM. 

54.  Ocean  boundaries  are  characterized  by  the  specification  of  a 
time  varying  water  suface  elevation  at  the  boundary.  Velocities  on  the 
ocean  boundary  are  then  computed  from  a  simplified  form  of  the  momentum 
equation  where  the  diffusive  terms  have  been  neglected.  One-sided 
differences  are  used  to  replace  derivatives  that  need  points  outside  the 
field.  As  an  example,  consider  the  computation  for  v  on  an  ocean 
boundary  that  lies  on  the  bottom  of  the  transformed  plan. 
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An  inspection  of  Equation  47  reveals  that  in  order  to  be  able  to  compute 
v“+1  ,  values  of  at  the  center  of  the  first  cell  must  be  known. 

These  are  determined  by  setting  them  equal  to  the  boundary  values  of  <f> 


(47) 
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I 

but  lagged  by  the  time  required  for  a  free  surface  gravity  wave  to 
traverse  the  distance  from  the  boundary  to  the  interior  point,  e.g., 
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55.  When  the  flow  is  directed  into  the  computational  field,  the  j 

,1 

boundary  condition  on  the  salinity  is  prescribed  as  that  of  the  ocean. 

| 

However,  when  the  flow  is  moving  out  of  the  computational  field,  the 
salinity  at  an  ocean  boundary  is  set  to  be  equal  to  its  value  at  the 
next  point  inside.  j 

56.  River  boundaries  are  characterized  by  the  specification  of 
the  velocity.  The  salinity  is  set  to  be  zero  and  the  water  surface 
elevation  at  the  center  of  a  river  boundary  cell  is  computed  as  in  any 
interior  cell. 
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PART  V:  MODEL  APPLICATIONS 

57.  In  order  to  demonstrate  the  versatility  of  VARM  in  its  ability 
to  model  flows  in  rather  general  multiply-connected  regions  containing 
both  river  and  ocean  boundaries,  three  applications  have  been  made  using 
the  physical  geometry  in  Figure  3. 

Generation  of  Boundary-Fitted  Coordinates 


58.  The  first  step  in  the  application  of  VAHM  is  the  generation 
of  the  boundary-fitted  coordinates.  This  is  accomplished  through  a 
coordinate  generation  code  developed  by  Thompson.  Output  from  the  co¬ 
ordinate  code  is  saved  on  a  file  for  subsequent  use  by  VAHM.  The  basic 
input  to  the  coordinate  code  is  the  specification  of  the  (x,y)  co¬ 
ordinates  of  the  boundary  points  noted  on  Figure  3.  Although  various 
degrees  of  coordinate  control  can  be  exercised,  the  boundary-fitted 
coordinates  shown  in  figure  3  were  computed  using  no  control.  Figure  4 
illustrates  the  actual  computational  grid  network  that  is  used  in  VAHM, 
where  velocities  are  computed  at  the  cell  corners  and  salinities  and 
water  elevations  at  the  cell  center.  However,  it  should  be  remembered 
that  VAHM  requires  that  the  (x,y)  coordinates  be  specified  at  not  only 
the  corners  and  center  of  a  computational  cell  but  also  on  the  cell 
faces.  The  reason  for  this  is  because  of  the  f.'ly  geometrically  con¬ 
servative  transformation  of  the  mass,  momentum  and  ,’alinity  equations  to 
be  solved.  With  such  a  transformation,  one  should  never  use  averaged 
values  of  the  geometrical  derivatives  since  this  can  result  in  the  loss 
of  conservation  of  the  properties  being  computed.  This  is  the  reason 
for  computing  the  coordinate  system  illustrated  in  Figure  3. 

59.  The  coordinate  system  plotted  in  Figure  3  was  the  third  at¬ 
tempt  at  generating  a  useful  grid  system.  Through  the  movement  of  bound¬ 
ary  points  and/or  coordinate  control  one  attempts  to  compute  boundary- 
fitted  coordinates  such  that  the  grid  spacing  does  not  vary  rapidly  and 
such  that  (£,n)  lines  never  approach  being  parallel  to  each  other.  The 
coordinate  system  in  Figure  3  satisfies  both  of  these  criteria  and  thus 
is  considered  to  be  adequate. 
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Flow  Through  Problem 


60.  Before  applying  the  model  for  the  case  of  time  varying  bound¬ 
ary  conditions,  various  "debugging"  applications  were  made.  Perhaps  the 
most  important  of  these  was  a  "flow  through"  problem.  In  a  flow  through 
test  all  velocities  are  set  to  be  equal,  but  non  zero,  (even  on  the 
walls),  the  water  surface  elevation  is  constant,  but  non  zero,  over  the 
field,  and  the  salinity  is  set  to  be  a  non  zero  constant  over  the  field. 
If  the  coding  is  correct  and  all  external  forces  have  been  set  to  zero, 
the  initial  state  should  never  change.  Such  tests  have  helped  to  correct 
many  errors  that  might  otherwise  have  gone  undetected. 

Case  1  -  Sloping  River 

61.  The  first  application  was  one  in  which  a  river  boundary  with 

a  constant  velocity  of  0.4  m/s  was  prescribed  at  the  top  with  the  water 

surface  elevation  at  the  bottom  being  held  constant  at  1.0  m.  The  bottom 

was  assumed  to  have  a  slope  of  0.005  m/grd  cell  and  the  initial  depth  was 

set  to  be  11.0  m.  The  initial  velocity  field  was  set  to  zero  as  was  the 

1/2 

salinity  concentration.  The  Chezy  coefficient  was  set  to  35  m  /s  and 
a  time  step  of  600  sec  was  prescribed.  These  plus  other  input  data  are 
presented  in  Table  1. 

62.  Three  separate  runs  were  made  in  which  the  influence  of  using 

a  form  of  Roache's  second  upwind  differencing  for  the  convective  terms 

in  the  momentum  equations  (C0NVEC  =  UPWIND)  as  opposed  to  centered 

differencing  (C0NVF.C  ~  CENTER)  and  the  influence  of  increasing  the  di- 

2  2 

agonal  components  of  the  eddy  viscosity  from  0.01  m  /s  to  10  m  s  were 
investigated . 

63.  Figures  5-8  illustrate  the  type  of  phenomena  that  can  occur 
when  using  centered  differences  to  represent  the  convective  terms.  After 
12  hours,  a  "zig  zag"  pattern  has  become  well  defined.  Figures  9-12 
demonstrate  the  effect  of  using  a  form  of  Roache's  second  upwind  differ¬ 
encing.  Although  a  slight  pattern  can  be  seen  after  12  hours,  it  isn't 
nearly  as  pronounced  as  when  using  centered  differences.  Figures  13-16 
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Table  1  -  Input  Data  to  VAHM 


Variable 

Case  1 

Case  2 

Case  3 

At  (sec) 

600 

600 

600 

D  (m2/s) 

XX 

0.01 

10 

10 

10 

D  (m2/s) 

yy 

0.01 

10 

10 

10 

D  (m2/s) 

xy 

0.0 

0.0 

0.0 

Dyx  (,n2/s) 

0.0 

0.0 

0.0 

Kx  (m2/s) 

0.01 

0.01 

0.01 

E  (m2/s) 

y 

0.01 

0.01 

0.01 

w 

CM 

w 

CJ 

35 

35 

35 

Initial  depth,  m 

11.0 

11.0 

11.0 

Bottom  slope/cell 

0.005 

0.0 

0.0 

CONVEC 

CENTER 

UPWIND 

UPWIND 

UPWIND 

Initial  Velocity,  m/s 

0.0 

0.0 

0.0 

Conv  Tolerance 

0.005 

0.005 

0.005 

i^itofctat^vt£»A‘.,,a  -*.....  ...-a... ..,. 


show  that  by  increasing  the  eddy  viscosity,  along  with  upwind  differenc¬ 
ing,  essentially  all  of  the  "roughness"  in  the  computed  velocity  field 
has  been  removed.  Figure  17  presents  a  time  history  of  the  water  surface 
profile  along  the  £  =  6  line. 


Case  2:  Closed  at  Tor 


Ocean  on  Bottom 


64.  The  second  application  was  for  the  case  of  a  closed  boundary 
at  the  top  and  an  ocean  boundary  on  the  bottom.  Once  again  the  initial 
velocity  and  salin  ty  fields  were  set  to  zero  and  the  initial  depth  was 
11.0  m.  Unlike  the  previous  application,  the  bottom  was  assumed  flat. 

65.  The  water  surface  elevation  curve,  relative  to  a  depth  of 

10  m,  presented  in  Figure  18  was  prescribed  at  the  ocean  boundary  along 

with  a  salt  concentration  of  30  ppt.  As  in  the  previous  application, 

CONVEC  =  UPWIND,  D  =  D  =10  m^/s  ,  At  =  600  sec  and  the  value  of 

xx  yy  l/2 

the  Chezy  coefficient  was  35  m  /s  (see  Table  1). 

66.  The  vector  plots  of  the  flow  field  presented  in  Figures  19-30 
illustrate  quite  clearly  the  effect  of  the  time  varying  ocean  boundary 
which  first  drives  water  into  the  field  with  water  flowing  out  on  the  ebb 
portion  of  the  tidal  cycle.  The  channelizing  effect  of  the  island  is 
also  quite  clearly  shown.  Figure  31  presents  a  time  history  of  the  water 
surface  profile  along  the  £  =  6  line,  whereas  Figures  32,  33  and  34 
are  plots  of  the  water  surface  at  particular  points. 

67.  From  an  inspection  of  Figures  26-30  it  can  be  seen  that  an 
oscillation  in  the  flow  field  has  developed  in  the  upper  portion  of  the 
modeled  area  when  the  flow  is  pushed  toward  the  boundary.  This  is 
probably  due  to  the  influence  of  the  upstream  boundary,  although  the 
tolerance  on  the  iterated  water  surface  may  also  be  a  factor.  The  con¬ 
vergence  tolerance  was  set  to  be  0.005  m  in  all  the  runs. 

68.  Figure  35  demonstrates  the  movement  of  the  salinity  field 
over  a  tidal  cycle.  As  the  flooding  cycle  of  the  tide  curve  is  experi¬ 
enced,  saline  water  at  a  concentration  of  30  ppt  moves  into  the  region. 

As  the  flow  reverses  at  the  ocean  boundary  (see  Figure  23),  the  salinity 
at  the  boundary  is  set  to  its  value  immediately  inside  to  reflect  an 
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outflow  boundary.  This  is  the  reason  for  the  decrease  in  salt  concentra¬ 
tion  at  the  boundary  after  5  hours.  Figure  36  gives  a  time  history  of 
the  salinity  at  a  point  about  12  miles  from  the  ocean  boundary. 


69.  The  third  application  was  identical  to  the  second  except  that 
a  river  boundary  with  a  constant  velocity  of  0.4  m/s  in  the  r)“co'11Ponent 
was  assumed  at  the  top  as  opposed  to  the  closed  boundary  in  Case  2. 

70.  Figures  37-49  present  "snap  shots"  of  the  computed  flow  field 
for  16  hours.  With  the  flow  field  initialized  to  zero  at  a  constant 
depth  of  11.0  m,  it  can  be  seen  that  the  influence  of  the  incoming  tide 
and  the  river  meet  after  about  4  hours.  As  in  the  previous  application, 
an  oscillatory  pattern  occurs  in  the  upper  portion  between  hours  8  and 
12.  However,  as  revealed  in  Figures  48  and  49,  this  irregularity  is  re¬ 
moved  as  the  influence  of  the  ebb  portion  of  the  tidal  cycle  is  felt  in 
the  upper  reach. 

71.  Figure  50  is  a  plot  of  the  time  history  of  the  water  surface 
profile  along  the  £  =  6  line  and  Figures  51-53  give  the  time  history 
of  elevations  at  particular  points.  Figure  54  presents  the  time  history 
of  the  salinity  at  a  point  about  12  miles  from  the  ocean  boundary.  Com¬ 
paring  Figure  54  with  Figure  36,  it  can  be  seen  that  the  salinities  are 
essentially  the  same  after  16  hours.  Therefore,  the  influence  of  the 
river  has  not  been  felt  in  the  lower  reach  after  this  length  of  time. 


Computing  Costs  and  Times 


72.  Previously  it  was  noted  that  the  major  advantage  of  the  FEM 
over  the  FDM  was  its  ability  to  more  accurately  handle  irregular  bound¬ 
aries,  whereas  its  major  disadvantages  were  increased  complexity  in 
coding  and  perhaps  increased  computational  costs.  A  finite  difference 
model  such  as  VAHM  which  makes  computations  on  a  boundary-fitted  co¬ 
ordinate  system  removes  the  boundary  representation  advantage  of  the  FEM 
in  many  problems.  Although  no  direct  comparison  can  be  made  between 
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VAHM  and  the  FEM  in  the  areas  of  complexity  of  coding  and  computational 
costs,  some  rather  general  comparisons  can  be  made. 

73.  Since  VAHM's  computational  scheme  retains  much  of  the  charac¬ 
ter  of  explicit  finite  difference  schemes,  it  would  appear  that  VAHM's 
coding  should  be  much  less  complicated  than  any  finite  element  model. 

With  simpler  coding,  future  modifications  should  be  much  easier  to  make, 
e.g.  flooding  of  boundaries,  higher  order  representation  of  the  advective 
terms,  etc. 

74.  An  approximate  comparison  of  computing  costs  can  be  made  with 

13 

a  vertically  averaged  finite  element  model  called  RMA-2  .  This  is  a 
flow  model  that  does  not  include  the  modeling  of  salinity  and  its  effect 
upon  the  flow  field.  The  model  was  developed  by  Resource  Management 
Associates  and  is  currently  being  used  by  the  Estuaries  Division  of  the 
Hydraulics  Laboratory  at  WES. 

75.  The  time  step  employed  in  the  previously  presented  runs  of 
VAHM  was  600  sec,  which  compares  with  a  time  step  of  perhaps  30  sec  that 
could  be  employed  in  a  fully  explicit  finite  difference  model.  With  a 
computational  grid  that  contains  363  velocity  points  (1365  coordinate 
points)  12  hours  of  computations  required  43  sec  of  CPU  time  at  a  cost  of 
$21  on  a  CRAY  I  computer  for  the  first  application  presented.  The  second 
and  third  applications  required  approximately  twice  as  much  CPU  time  at 
about  twice  the  cost.  However,  later  experimentation  with  VAHM  revealed 
that  stable  computations  could  be  achieved  using  a  time  step  of  1500  sec. 
Therefore,  if  such  a  time  step  had  been  used  in  the  cases  presented  the 
costs  would  have  been  reduced  by  a  factor  of  about  2.5.  The  increased 
time  of  the  last  two  applications  was  because  of  the  more  rapidly  varying 
water  surface.  Only  one  or  two  iterations  each  time  step  were  required 
to  achieve  a  convergence  tolerance  of  0.005  m  in  the  first  case,  whereas 
an  average  of  seven  or  eight  iterations  were  required  in  the  last  two 
applications . 

76.  As  a  comparison,  RMA-2  applications  to  grids  containing  ap¬ 
proximately  the  same  number  of  net  points  for  a  12  hour  tidal  cycle,  using 

14 

an  1800  sec  time  step,  cost  about  $40  on  the  same  CRAY  I  computer. 

However,  it  should  be  remembered  that  RMA-2  makes  computations  for  only 
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the  flow  field,  whereas  VAHM  also  computes  the  salinity  field  and  its 
coupling  with  the  flow  field  through  a  relationship  with  the  water 
density.  Based  upon  these  approximate  costs,  it  woul'i  appear  that  the 
FEM,  as  reflected  by  RMA-2  costs,  is  about  3  times  more  expensive  than 
VAHM  for  tidal  problems  and  perhaps  6  times  more  expensive  for  river 
problems.  As  a  final  note,  it  is  believed  that  after  VAHM  has  been 
"cleaned  up"  to  better  utilize  the  vector  processing  features  of  the 
CRAY  I  its  computational  costs  will  decrease  significantly. 


PART  VI:  SUMMARY  AND  RECOMMENDATIONS 


77.  A  numerical  model  for  computing  vertically  averaged  velocities 
and  salinity  plus  water  surface  elevations  has  been  developed.  By  em¬ 
ploying  the  concept  of  boundary-fitted  coordinates,  irregular  boundaries 
can  be  accurately  modeled  in  either  simply  or  multiply-connected  regions. 
Even  though  the  numerical  grid  is  a  nonorthogonal  curvilinear  grid  in  the 
physical  region  being  modeled,  all  numerical  computations  are  carried  out 
in  a  transformed  rectangular  grid  with  square  grid  spacing. 

78.  A  feature  of  the  model  is  the  particular  solution  technique 
employed  to  numerically  solve  the  governing  equations.  A  combination 
implicit-explicit  finite  difference  scheme,  patterned  after  work  by 
Edinger  and  Buchak^  in  their  development  of  a  laterally  averaged  reser¬ 
voir  hydrodynamic  model,  has  been  developed  to  remove  the  speed  of  a 
gravity  wave  from  stability  restrictions  on  the  computational  time  step 
while  still  retaining  some  of  the  advantages  of  explicit  schemes.  With 
such  a  scheme,  the  water  surface  elevation  is  computed  implicitly  using 
the  Accelerated  Gauss-Seidel  solution  technique  while  the  velocities  and 
salinity  are  computed  in  an  explicit  fashion. 

79.  The  model  has  been  developed  for  general  applications.  Any 
number  of  river  and/or  ocean  boundaries  can  be  arbitrarily  located  on 
the  transformed  rectangular  plane,  as  can  the  placement  of  islands  in 
the  interior  of  the  computational  field.  Even  though  a  great  deal  of 
generality  does  exist,  there  are  restrictions.  For  example,  only  no-slip 
boundary  conditions  are  currently  treated  at  solid  boundaries  and  no 
flooding  of  those  boundaries  is  allowed. 

80.  Although  VAHM  has  been  developed  to  the  point  where  results 
from  the  test  applications  presented  are  encouraging,  additional  work  is 
needed  before  VAHM  can  be  considered  fully  operational.  Recommendations 
for  additional  development  are  listed  below. 

-  In  order  to  expand  VAHM’s  capabilities  into  the  water  quality 
area,  it  is  necessary  to  devise  a  scheme  for  solving  the 
transport  equation  that  accurately  transports  a  "spike"  concen¬ 
tration  distribution.  The  present  scheme  employed  in  VAHM  for 
computing  salinity  is  sufficient  when  distributions  are  fairly 
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smooth.  However,  it  will  not  do  a  good  job  on  a  spike  distri¬ 
bution.  Therefore,  a  major  task  to  be  accomplished  is  the 
modification  of  VAHM  to  allow  such  computations  to  be  ac¬ 
curately  made. 

-  As  previously  discussed,  at  the  present  time  only  no-slip 
conditions  are  allowed  at  solid  boundaries.  Unless  a  small 
grid  spacing  is  used  near  the  boundaries,  slip  conditions  may 
be  more  appropriate.  The  slip  boundary  conditions  will  be 
determined  ay  requiring  the  normal  component  of  the  velocity 
and  the  vorticity  to  be  zero  at.  a  wall.  Many  of  the  checks 
in  the  basic  model  have  been  coded  to  allow  for  slip  condi¬ 
tions;  however,  the  slip  subroutine  remains  to  be  developed. 

-  The  capability  of  handling  flooding  boundaries  is  needed 

in  VAHM.  Some  ideas  for  incorporating  such  a  capability  into 
VAHM  have  been  considered  in  the  basic  coding. 

-  VAHM  uses  the  Accelerated  Gauss-Seidel  solution  technique  to 
implicitly  compute  the  water  surface  elevation.  At  the  pres¬ 
ent  time,  a  constant  acceleration  parameter  is  employed.  The 
use  of  variable  acceleration  parameters  for  the  purpose  of 
speeding  up  the  computations  should  be  investigated. 

-  At  the  present  time,  a  2D  vector  plotting  program  developed 
by  S.  A.  Adamec  of  the  Hydraulics  Laboratory  at  WES  has  been 
coupled  with  the  grid  generation  code  and  VAHM  to  provide 
plots  of  the  velocity  field.  Additional  plotting  capability 
needs  to  be  coupled  with  VAHM. 

-  As  noted,  the  basic  coding  has  been  written  to  allow  for  an 
extremely  general  representation  of  a  physical  problem,  e.g, 
specification  of  islands,  r.’.ver  inlets,  etc.  However,  many 
of  these  options  have  not  been  "debugged."  In  addition, 
although  VAHM  is  being  run  on  a  CRAY-I  computer  no  attempt 
at  "cleaning  up"  the  code  to  take  advantage  of  the  CRAY's 
vector  processing  has  been  made. 
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Figure  3.  Boundary  fitted  coordinate  system  for  a  multiple 
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Figure  4.  Computational  grid  for  computing  velocities 
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Figure  5.  Velocity  field  after  3  hours,  CONVEC  »  CENTER 
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Figure  6.  Velocity  field  after  6  hours,  CONVEC  -  CENTER 
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Figure  7.  Velocity  field  after  9  hours,  CONVEC  «  CENTER 
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Figure  8.  Velocity  field  after  12  hours,  CONVEC  -  CENTER 
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Figure  10.  Velocity  field  after  6  hours,  CONVEC  ■  UPWIND 
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Figure  11.  Velocity  field  after  9  hours,  CONVEC  -  UPWIND 
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Figure  12.  Velocity  field  after  12  hours,  CONVEC  -  UPWIND 


turn 

4 
" . . 

*4****** 

4*4444*4^ 

**4444444 

♦4*44444* 

44*4*4444 

4*44*4*44 

***4 

♦»  **','$7 

»  r»,  4 i 7 

'*  *  4  *  /  * 

*  •  *  *  «  * 


VELOCITY 

WCCTOA 


EXCEEDS  PLOT  LIMIT 


3.2:  KM/ IN 
9.66  KM/ IN 


Figure  14.  Velocity  field  after  6  hours,  CONVEC  ■  UPWIND,  • 
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Figure  15.  Velocity  field  after  9  hours,  CONVEC  -  UPWIND,  D 


Figure  17.  Water  surface  profile  along  £  «  6  line 
for  a  sloping  river  with  control  at  downstresun 
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Figure  20.  Velocity  field  after  2  hours  with  an  ocean  boundary 
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Figure  21.  Velocity  field  after  3  hours  with  an  ocean  boundary 
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Figure  22.  Velocity  field  after  A  hours  with  an  ocean  boundary 
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Figure  23.  Velocity  field  after 
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Figure  24.  Velocity  field  after  6  hours  with  an  ocean  boundary 
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Figure  25.  Velocity  field  after  7  hours  with  an  ocean  boundary 
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Figure  26.  Velocity  field  after  8  hours  with  an  ocean  boundary 
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Figure  28.  Velocity  field  after  10  hours  with  an  ocean  boundary 
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Figure  29.  Velocity  field  after  11  hours  with  an  ocean  boundary 
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Figure  30.  Velocity  field  after  12  hours  with  an  ocean  boundary 
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Figure  36.  Salinity  at  £  =  6,  1 
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TO 

,•%***» 

1.8 

VELOCITY 

VECTOR 

SCALE 

1.1 

_  m- 

cup* ) 

9T  LIMIT 

«t*,f'**  , 

iff'"* 

,*• 

xs  ■ 

¥8  • 

UKtn  ru 

3.22  KR/IH 

s.M  km/in 

v‘v 

'  y » 


4 


nnnit’t 

iititmt 


f  f  t  ♦  t  ♦  t  t  i 


Figure  37.  Velocity  field  after  1  hour  with  an  ocean  and  a  river  boundary 
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Figure  38.  Velocity  field  after  2  hours  with  an  ocean  and  a  river  boundary 
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Figure  39.  Velocity  field  after  3  hours  with  an  ocean  and  a  river  boundary 


'3  ’■  tA- '  V  J  tf  I*  If- \  i>2- !•  . T  f'IVi «' }0  ■•’Vi* fii  I1*  h  Aif  W  ‘  1 


UCLOC I TV 

UCCTOK 

SCALE 

- - —  * **  tN*S> 

l.t 

<+——  EXCEEDS  NIOT  LINIT 

MS  •  3.22  in/ IN 
VS  •  9.SS  IN/ IN 


n»«m» 

Mf«V 


W;:  ttt 


-it 

I  I  7  ]  f  ■  a  |  u 

i  ii !  m  n 


Figure  40.  Velocity  field  after  4  hours  with  an  ocean  and  a  river  boundary 
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Figure  41.  Velocity  field  after  5  hours  with  an  ocean  and  a  river  boundary 


Figure  42.  Velocity  field  after  6  hours  with  an  ocean  and  a  river  boundary 
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Figure  43.  Velocity  field  after  7  hours  with  an  ocean  and  a  river  boundary 
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Figure  44.  Velocity  field  after  8  hours  with  an  ocean  and  a  river  boundary 
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Figure  45.  Velocity  field  after  9  hours  with  an  ocean  and  a  river  boundary 
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Figure  46.  Velocity  field  after  10  hours  with  an  ocean  and  a 

river  boundary 
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Figure  47.  Velocity  field  after  12  hours  with  an  ocean  and 

river  boundary 
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Figure  48.  Velocity  field  after  14  hours  with  an  ocean  and  a 

river  boundary 
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Figure  49.  Velocity  field  after  16  hours  with  an  ocean  and 

river  boundary 
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APPENDIX  A:  NOTATION 


A.  . 

ij 

A.  . 
ij 

a 


D 

D 


,D 

xx  yy 
,D 

xy  yx 

E  ,E 
x  y 


f  f  ft,f 

x  y  V  n 


g 

h 

j 


p.Q 

p 


T 

At 


u,v,w 


u.  ,u . ,u. 
i  J  k 


u . 

l 

u! 

1 

n. 

i 


u! 

i 

v 

w 

x,y,z 

4,1 

Ax, Ay, A£, An 


Eddy  diffusivity  tensor 
Eddy  dispersion  tensor 
Constant 

Chezy  coefficient 
Constant,  distance 
Molecular  diffusivity 

Diagonal  components  of  eddy  viscosity  tensor 
Off  diagonal  components  of  eddy  viscosity  tensor 
Components  of  eddy  dispersion  tensor 
Arbitrary  function 
Derivatives 

Acceleration  of  gravity 
Water  depth 

Jacobian  of  the  transformation 
Coordinate  control  functions 
Pressure 

Atmospheric  pressure 

Cylindrical  coordinates 

Salinity 

Temperature 

Time  step 

Components  of  velocity 
Tensor  notation  for  velocity 

Time  averaged  velocity 

Random  time  varying  component  of  velocity 
Time  and  depth  averaged  velocity 

Random  depth  varying  component  of  time  averaged  velocity 
Wind  speed 

Cartesian  coordinates 
Boundary-fitted  coordinates 
Spatial  grid  steps 


A1 


3/9x. ,3/3x. 
i  J 

X  ,T 
s  s 
x  y 

tb  *tb 

x  y 


Water  density 

Reference  water  density 

Density  of  air 

Molecular  viscosity 

Turbulent  viscosity  tensor 

Eddy  dispersion  viscosity  tensor 

Water  surface  elevation;  arbitrary  variable 

Stress  tensor 

Cylic  tensor 

Kronecker  delta 

Time  derivative 

Space  derivatives 

Components  of  bottom  shear  stress^^ 
Components  of  bottom  shear  stress Jq 
Wind  direction 

Latitude  of  center  of  modeled  area 
Earth's  angular  velocity 
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In  accordance  with  letter  from  DAEN-RDC,  DaEN-ASI  dated 
22  July  1977,  Subject:  Facsimile  Catalog  Cards  for 
Laboratory  Technical  Publications,  a  facsimile  catalog 
card  in  Library  of  Congress  MARC  format  is  reproduced 
below. 


Johnson,  Billy  H 

VAHM  -  A  vertically  averaged  hydrodynamic  model 
using  boundary-fitted  coordinates  /  by  Billy  H.  Johnson. 
Vicksburg,  Miss.  :  U.  S.  Waterways  Experiment  Station  ; 
Springfield,  Va.  :  available  from  National  Technical 
Information  Service,  1980. 

52,  [56]  p.  :  ill.  ;  27  cm.  (Miscellaneous  paper  - 
U.  S.  Army  Engineer  Waterways  Experiment  Station  ;  HL-80-3) 
Prepared  for  Assistant  Secretary  of  the  Army  (R&D), 
Department  of  the  Army,  Washington,  D.  C.,  under  Project 
UAO0IIOIA9ID. 

References:  p.  51-52. 

1.  Computerized  models.  2.  Coordinates.  3.  Hydrodynamics, 

4.  Mathematical  models.  5.  Numerical  analysis,  6.  Salinity. 
7.  VAHM  (Vertically  Averaged  Hydrodynamic  Model).  I.  United 
States.  Assistant  Secretary  of  the  Army  (Research  and 
Development).  II.  Series:  United  States.  Waterways  Experiment 
Station,  Vicksburg,  Miss.  Miscellaneous  paper  •,  HL-80-3. 
TA7.W3**m  no.  HL-80-3 


